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13.  ABSTRACT  (Maximum  200  words) 


This  project  aimed  first  at  obtaining  a  better  understanding  of  numericai  methods  for  hyperbolic  conservation  laws.  Several 
significant  advances  were  obtained.  In  particular,  the  interplay  between  the  size,  uniformity  and  topology  of  numerical  gnds  and  the 
a^racy  of  numerical  solutions  based  on  those  grids  has  been  uncovered.  This  type  of  results  in  unique  for  problems  with  low 
regularity.  Those  a  priori  estimates  have  led  to  the  first  explanation  of  supraconvergence  for  conservation  laws.  The  possrbritty  of 
using  earlier  results,  a  posteriori  estimates,  for  the  design  of  adaptive  methods  with  full  mathematical  backing  was  also  investigated. 


In  a  second  part  numerical  methods  for  related  problems  were  considered.  For  systems  of  Hamilton-Jacobi  equations,  the  notion  of 
viscosity  solutiori  cannot  be  applied.  In  spite  of  this,  competitive  numerical  methods  were  constructed  and  successfully  used.  The 
flow  of  granular  materials,  a  highly  important  but  ill-understood  problem,  was  also  considered.  Those  problems  have  several 
nonstandard  features.  Computational  results  based  on  the  first  implementation,  in  that  field,  of  higher  order  numencal  (Discontinuous 
Galerkin)  methods  were  obtained.  A  code  is  currently  being  tested  by  engineers  in  an  industrial  setting. 


1.  LIST  OF  FIGURES: 


(1)  Figure  1:  Geometry  of  the  hopper  problem,  p.5. 

(2)  Figure  2:  Structure  of  the  stress  field:  rarefaction  case,  p.6. 

(3)  Figure  3:  Structure  of  the  stress  field:  shock  case,  p.6. 

(4)  Figure  4:  Graphic  user  interface  for  the  powder  problem,  p.8. 


2.  STATEMENT  OF  THE  PROBLEMS  STUDIED: 

The  various  goals  of  the  project  were  clearly  stated  in  the  main  proposal.  Seven  of  them,  widely  different 
in  length  and  difficulties  were  explicitly  listed;  they  were 

(1)  derivation  of  optimal  error  estimates  for  multidimensional  monotone  schemes  for  conservation  laws, 

(2)  generalization  to  schemes  that  do  not  possess  splitting  fluxes, 

(3)  generalization  to  higher  order  methods, 

(4)  study  of  the  possibility,  or  lack  thereof,  of  generalization  to  systems  of  conservation  laws, 

(5)  derivation  of  general  a  posteriori  estimates, 

(6)  design  and  implementation  of  mathematically  sound  adaptive  strategies  for  hyperbolic  problems, 

(7)  applications  to  granular  flows  and  other  problems  (in  particular  Hamilton- Jacobi  equations,  even 
though  those  were  not  explicitly  mentioned  in  the  proposal). 


3.  SUMMARY  OF  THE  MOST  IMPORTANT  RESULTS: 

Significant  advances  were  obtained  for  most,  but  not  all  of  the  above  problems.  People  who  collaborated 
with  the  PI  on  questions  covered  by  the  project  at  one  point  or  another,  but  were  not  directly  supported  by 
it,  include  Profs.  G.  Amiez,  (U.  of  Franche-Comte,  France),  V.  Alexiades  (U.  of  Tennessee),  B.  Cockburn  (U. 
of  Minnesota),  D.  Schaeffer  (Duke  University)  and  M.  Shearer  (North  Carolina  State  University),  as  well  as 
K.S.  Lowder  (former  graduate  student,  Gremaud  adviser,  U.S.  Airforce  Academy),  N.R.  Ide  (former  graduate 
student,  Gremaud  adviser,  Lincoln  Lab,  MIT),  J.V.  Matthews  (current  Ph.D.  student,  Gremaud  adviser), 
and  W.  Turner  (former  graduate  student,  Gremaud  adviser,  NCSU).  We  have  also  worked  in  consultation 
with  Engineers  at  Jenike  k  Johanson,  Inc.,  Westford,  MA.,  and  in  particular  T.A.  Royal,  Vice  President,  to 
develop  efficient  and  robust  solvers  for  various  types  of  problems  related  to  granular  flows. 

The  accomplishments  corresponding  to  the  above  categories  are  now  reviewed.  In  a  first  part,  the  work 
pertaining  to  error  estimates  for  conservation  laws  is  discussed  (parts  (1)  through  (6)).  A  second  part  deals 
with  the  advances  related  to  both  Hamilton- Jacobi  equations  and  granular  materials.  Finally,  some  remarks 
on  ongoing  and  future  work  are  offered. 


Error  estimates  for  numerical  conservation  laws  (parts  (l)-(6)). 

The  goal  was  here  to  estimate  the  error  between  the  weak  entropy  solution  u  to  a  generic  scalar  conser¬ 


vation  law 


-h  V  •  f(v)  =  0, 


t'(-,0)  =  t;o. 

and  a  numerical  solution  u  constructed  from  as  general  a  method  as  possible.  Several  aspects  of  those  types 
of  hyperbolic  problems  makes  this  a  challenging  task,  namely  the  presence  of  shocks  and  the  existence  of 
many  weak  solutions,  among  which  only  the  entropy  solution  is  believed  to  have  a  physical  meaning. 

Our  main  contribution  is  the  derivation  of  a  priori  estimates,  i.e.,  estimates  written  only  in  terms  of  the 
exact  solution.  Those  estimates  are  of  the  form 


e{T)<^v,T), 

where  e(t)  denotes  the  L^-error  at  time  t.  Our  results  are  the  first  and  only  rigorous  error  estimates  of  this 
type  for  nonlinear  conservation  laws.  All  the  other  results  we  are  aware  of  are  of  the  form 

e(T)  <  ^(v,u,T), 


Typeset  by  ^,^45-T]eX 


i.e.,  they  also  involve  the  approximate  solution  u.  If  does  not  involve  v  but  only  u,  the  latter  estimate  is 
referred  to  as  an  a  posteriori  error  estimate. 

At  the  beginning  of  this  proposal,  virtually  every  error  estimate  was  based  on  an  original  approach  by 
Kuznetsov  (1976).  In  a  joint  work  with  B.  Cockburn  (1996)  that  took  place  slightly  before  the  start  of 
the  present  project,  we  identified  Kuztnetsov’s  point  of  view  as  an  a  posteriori  approach.  This  was  then 
applied  to  the  derivation  of  a  posteriori  error  estimates  for  both  shock  capturing  streamline  diffusion  methods 
and  shock  capturing  discontinuous  Galerkin  methods.  A  posteriori  estimates  are  very  useful  for  practical 
computations  if  they  are  properly  used  as  part  of  an  adaptive  strategy.  Part  of  the  present  project,  see  part 
(6)  above,  dealt  with  those  very  issues.  Essentially,  our  a  posteriori  estimates  are  as  follows 

e(T)  <  e(0)  +  C\\  dtHniu)  +  V  •  /(t,  a:,  nA(u))  ||i.(o,T,L>) 

where  C  is  a  numerical  constant,  11/,  an  interpolation  operator  into  a  finite  dimensional  space  of  continuous, 
and  thus  5tn/,(M)  +  V  •  f{t,  x, !!/,(«))  is  the  residual.  The  feasibility  of  using  such  a  result  as  the  basis  for 
an  adaptive  strategy  was  studied  for  instance  in  K.  Lowder,  Mesh  Optimization  for  Conservation  Laws, 
Master’s  Project,  Dept,  of  Mathematics,  NCSU,  Gremaud  Adviser  (1996).  The  point  of  view  taken  there 
was  to  obtain  a  “close  to  optimal”  mesh  at  each  time,  based  on  the  error  estimate.  The  results  were  mixed. 
In  the  simplest  cases  tried,  the  grids  were  successfully  adapted  to  the  features  of  the  solution.  However, 
first,  it  came  at  considerable  computational  cost.  Second,  for  more  involved  cases,  the  construction  of  the 
residual,  for  instance  the  choice  of  a  proper  interpolation  operator  H/,  in  the  above  expression,  played  a 
crucial  role.  We  were  not  able  to  find  a  general  way  of  constructing  such  a  operator  that  would  work  for 
“all”  cases.  Here  again,  having  to  deal  with  possibly  discontinuous  solutions  is  a  significant  complicating 
factor.  The  algorithms  obtained  were  not  competitive  with  most  of  the  common  heuristic  methods  in  this 
field,  as  far  as  adaptivity  is  concerned. 

A  natural  next  step  consisted  then  in  improving  our  understanding  of  the  influence  of  the  grid  on  the 
accuracy  of  the  numerical  solution.  Indeed,  an  adaptive  strategy  is  likely  to  change  the  size,  the  uniformity 
and  the  topology  of  the  grid.  The  best  way  of  studying  those  effects  is  through  a  priori  error  estimates.  This 
project  has  led  to  very  significant  improvements  in  that  respect. 

It  has  been  known  from  quite  some  time  that  if  the  numerical  fluxes  do  not  properly  take  into  account  the 
irregularities  of  the  grid,  a  loss  of  consistency  occurs.  This  can  be  shown  on  very  simple  examples,  for  which 
the  truncation  error  does  not  go  to  zero,  see  e.g.  [8]  below.  The  fact  is,  that  in  most  practical  codes,  numerical 
fluxes  are  not  modified  to  take  into  account  local  mesh  variations.  Yet,  those  codes  perform  quite  well,  and 
not  only  converge  to  the  correct  solution,  but  may  in  fact  do  so  without  loss  in  the  order  of  convergence,  this, 
in  spite  of  the  loss  of  consistency.  This  supraconvergence  phenomenon  remained  unexplained  for  nonlinear 
conservation  laws  until  the  publication  of  [8],  and  [7]  after  that.  The  approach  is  based  on  a  modification  of 
Kuznetsov’s  original  approximation  theory.  The  point  of  view  taken  is  essentially  the  dual  of  Kuznetsov’s, 
and  does  lead  to  a  sound  a  priori  theory.  This  was  applied  in  [8]  to  monotone  schemes  on  nonuniform 
Cartesian  grids,  and  in  [8],  to  monotone  on  general  multidimensional  grids.  An  interesting  point  is  related 
to  the  importance  of  the  topology  of  the  grid.  It  is  known  (Sanders,  1983)  that  for  monotone  schemes  on 
Cartesian  grids,  uniform  or  not,  the  total  variation  of  the  numerical  solution  does  not  increase  in  time. 
This  is  not  true  for  non-Cartesian  grids,  even  uniform  ones.  Indeed,  in  [7],  a  simple  example  is  presented 
on  a  uniform  grid  consisting  of  equilateral  triangles,  and  hence  non-Cartesian,  for  which  even  the  Lax- 
Friedrichs  scheme  presents  an  increase  in  total  variation.  This  shows  that  the  case  of  non-Cartesian  grids 
is  intrinsically,  not  just  technically,  more  complex.  The  estimates  in  [7]  are  given  in  terms  of  quantities 
measuring,  for  general  multidimensional  grids,  the  departure  from  Cartesianity  and  uniformity.  This  work 
is  the  first  rigorous  explanation  of  the  abovementioned  supraconvergence  phenomenon.  In  essence,  what 
happens  is  that  the  possible  lack  of  consistency  of  the  numerical  methods  is  in  fact  compensated  providing 
two  fundamental  properties  are  satisfied:  first,  the  methods  are  conservative,  second,  the  numerical  fluxes 
are  consistent  (not  to  be  confused  with  the  consistency  of  the  entire  method). 
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The  above  approach  can  in  principle  be  applied  to  higher  order  methods.  It  is  in  fact  the  only  theoretical 
approach  we  know  for  which  this  statement  can  be  made  in  a  reasonably  general  sense.  However,  the  shear 
technical  complexity  is  enormous.  This  point  remains  under  study.  Another  aspect  was  the  possibility,  or 
lack  thereof,  of  a  generalization  of  the  above  results  to  systems  of  hyperbolic  conservation  laws.  Here,  the 
difficulties  are  not  technical,  but  conceptual,  since  in  fact  the  whole  notion  of  entropy  solution  becomes 
problematic  for  systems.  The  Numerical  Analysis  of  such  problems  remains  essentially  an  open  problem. 
It  would  appear  that  what  is  needed  is  not  some  kind  of  incremental  improvement,  built  on  the  available 
theory,  but  a  complete  breakthrough.  There  have  been,  however,  progresses  made  during  the  period  of  this 
project  related  to  the  understanding  of  systems  of  first  order  nonlinear  Partial  Differential  Equations,  see  in 
particular  the  work  below  in  systems  of  Hamilton- Jacobi  equations. 

Applications:  Hamilton- Jacobi  equations  and  granular  materials. 

As  mentioned  above,  the  case  of  systems  of  conservation  laws  is  still  largely  an  open  problem.  Recent 
results  have  however  been  obtained  for  Hamilton- Jacobi  equations,  both  in  the  scalar  and  vectorial  cases 
(Dacorogna-Marcellini,  1996).  Essentially,  those  results  apply  to  cases  for  which  the  usual  viscosity  approach 
is  either  ill-suited  or  not  applicable.  An  important  application  is  solid-solid  phase  transition  problems,  which 
are  usually  expressed  as  multi- well  potential  problems.  Those  can  be  recast  as  systems  of  Hamilton- Jacobi 
equations,  an  entirely  new  approach. 

We  have  proposed  and  tested  a  numerical  approach  for  the  calculation  of  solutions  to  Hamilton- Jacobi 
equations  under  a  viscoelasticity/capillarity  criterion.  The  (very  weak)  notion  of  “satisfying  the  boundary 
conditions  in  the  viscous  sense”  as  well  as  the  nonexistence  of  the  very  notion  of  viscosity  solution  in  the 
vectorial  case  prevented  us  from  using  the  classical  viscosity  criterion.  The  challenge  was  to  construct 
methods  that  would  convergence  to  a  properly  selected  solution  in  the  scalar  and  in  the  vectorial  case.  The 
idea  of  the  methods  consists  in  discretizing  the  original  PDEs  by  finite  difference  formulas  of  sufficiently 
high  order  as  to  guarantee  the  model  equations  to  be  identical  to  the  PDEs.  Anything  short  of  that  would 
lead,  and  indeed  did  lead,  to  an  effective  modification  of  the  selection  criterion.  The  approach  is  used  in 
conjunction  with  a  pseudotime  integrator  of  the  BDF  type.  Before  this  work,  almost  all  the  results  dealing 
with  the  use  of  the  above  selection  criterion  were  for  scalar  equations  in  the  Id  case.  We  have  not  only 
treated  multidimensional  scalar  problems,  but  also  multidimensional  vectorial  ones.  We  are  not  aware  of 
other  numerical  results  about  vectorial  Hamilton- Jacobi  equations.  Practically  speaking,  the  advantage 
of  the  method  is  as  follows.  For  typical  solid-solid  phase  change  problems,  the  traditional  approach  is  to 
construct  an  energy  functional  (always  nonconvex,  unless  the  problem  is  trivial),  discretize  the  corresponding 
variational  problem,  and  then  solve  somehow  the  corresponding  discrete  strongly  nonconvex  minimization 
problem.  As  is  well  known,  nonconvex  minimization  is  hard,  and  the  methods  likely  to  be  used,  at  least  in 
the  case  of  solid-solid  phase  change  problems,  are  likely  to  be  ad  hoc,  and  void  of  any  physical  justification. 
The  aim  here  was  to  bypass  the  minimization  problem  by  solving  to  steady  state  a  physically  meaningful 
problem  (time  dependent  Hamilton- Jacobi  system).  The  feasibility  of  the  approach  was  illustrated  in  [3]. 

The  last  part  of  the  project  involves  the  study  of  an  important  type  of  application:  the  study  of  the  flow 
of  granular  materials.  In  the  general  subject  of  industrial  flows,  it  can  be  argued  that  granular  flows  are  as 
important  as  fluid/gas  flows  (see  agriculture,  chemical,  construction  and  pharmaceutical  industries  to  name 
but  a  few).  In  spite  of  this,  those  problems  have  received  very  little  attention.  Part  of  this  project  deals 
with  those  problems  (see  part  4,  project  abstract,  p.  C-1.). 

One  of  the  first  and  main  problems  one  encounters  when  dealing  with  granular  materials  is  to  obtain  a 
proper  mathematical  model.  Several  models  exist,  that  appear  to  capture  well  some,  but  in  general  not  all, 
aspects  of  the  phenomena  under  consideration.  Our  approach  is  twofold.  First,  we  aim  at  obtaining,  through 
careful  Analysis  and  Numerical  Analysis  of  existing  models,  efficient  and  reliable  computational  tools  that 
will  allow  for  meaningful  investigations  the  models  themselves,  through  comparisons  with  experiments  for 
instance.  Second,  and  at  a  less  fundamental  level,  having  such  computational  tools  represents  a  breakthrough 
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in  this  field,  and  open  new  opportunities  in  terms  of  numerical  simulations  in  industrial  settings. 

So  far,  our  study  has  been  restricted  to  established  hopper  flows.  This  constitutes  in  our  view  a  first, 
meaningful,  step  toward  the  study  of  more  general  problems,  as  discussed  below.  Equations  governing  the 
time  dependent  flow  of  granular  material  under  gravity  can  be  derived  and  analyzed.  However,  those  are 
found  to  be  linearly  illposed  in  most  cases  of  practical  interest.  To  our  knowledge,  the  situation  is  not 
fully  understood,  mathematically  or  otherwise.  In  practice,  strongly  time-dependent  problems  are  usually 
observed  in  conjunction  with  funnel  Hows,  i.e.,  flows  for  which  the  motion  is  essentially  restricted  to  the 
central  part  of  the  silo.  Our  study  has  dealt,  so  far,  exclusively  with  mass  Hows,  i.e.,  flows  for  which  all  the 
material  is  mobilized.  In  this  context,  established,  steady  state  flows  can  be  observed. 

In  a  nutshell,  our  contributions  have  been 

(1)  construction  and  implementation  of  cellular  automaton  model  for  hopper  flows, 

(2)  investigation  of  the  properties  and  stability  of  radial  stress  fields  for  various  plasticity  models, 

(3)  investigation  of  extensions  of  the  notion  of  radial  solutions,  and  their  possible  use  as  design  tools  for 
flow  corrective  inserts, 

(4)  first  implementation  of  a  highly  accurate  method  (Discontinuous  Galerkin  method)  for  the  compu¬ 
tation  of  hopper  flows, 

(5)  computation  of  “switch  stresses”  due  to  changes  of  the  hopper’s  wall  angle. 

We  now  proceed  to  describe  in  more  details  the  above  points. 

Cellular  automaton.  A  simple  cellular  automaton  model  was  designed.  The  idea  is  essentially  to  work 
with  discrete  particles,  that  can  only  be  in  “states”  taking  values  in  some  discrete  spaces.  In  spite  of  its 
simplicity,  this  type  of  models  allowed  us  to  predict  structure  formations  reminiscent  of  the  density  waves 
observed  experimentally  with  some  hopper  flows,  see  [9]  for  details.  However,  the  approach  also  suffers  from 
serious  drawbacks.  In  particular,  by  its  very  nature,  all  the  information  available  from  such  simulations  is 
at  the  particle  level.  It  is  not  always  obvious  and/or  convenient  to  translate  this  into  quantities  that  can 
be  experimentally  measured.  Further,  an  obvious  bottleneck  is  the  number  of  particles  that  can  be  treated. 
The  continuum  approach,  that  we  have  adopted  since  then,  has  none  of  those  disadvantages. 

Radial  solutions.  In  some  cases  of  industrial  importance,  such  as  axisymmetric  containers  (silos), 
similarity  solutions  (radial  solutions)  can  be  constructed.  Most  of  the  existing  work  in  this  field  is  in  fact 
based  the  use  of  those  solutions.  This  is  especially  true  of  the  Engineering  literature.  In  [6],  those  solutions 
were  constructed  for  two  plasticity  models,  namely  for  Mohr-Coulomb  and  von  Mises  yield  conditions. 
Comparisons  of  the  stress  and  velocity  fields  for  the  two  models  gave  a  better  picture  of  their  respective 
properties.  Further,  in  the  simplest  of  the  two  cases  (Mohr-Coulomb),  a  numerical  stability  study  was 
performed.  Thanks  to  the  use  of  a  highly  accurate  method  (pseudospectral  Chebyshev  collocation  method), 
a  fairly  complete  picture  of  the  influence  of  the  material  parameters,  such  as  internal  friction  and  wall  friction, 
on  the  stability  was  uncovered. 

Extension  of  the  notion  of  radial  solutions.  The  notion  of  radial  solution  has  several  obvious 
shortcomings.  Prime  among  those  is  the  fact  that  it  only  exists  in  radially  symmetric  domains.  As  an 
attempt  to  alleviate  this  limitation,  we  have  stndied  the  possibility  of  considering  the  radial  solution  as  the 
first  term  in  an  expansion  of  the  following  type 

r(r,  0)  =  rT\e)  +  r^T^{9)  +  ... 

v(r,  0)  =  — u®(0)  -h  v^(0)  -h  . . . . 
r 

where  T  and  v  stand  for  generic  stress  and  velocity  components  respectively.  The  zero-th  order  terms 
correspond  to  the  radial  solution.  This  approach  was  explored  in  [6]  were  it  was  applied  to  non  radially 
symmetric  domains,  and  specifically  to  the  determination  of  flow  corrective  inserts.  Although  the  approach 
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was  successful,  it  also  brought  to  light  several  limitations  of  even  this  type  of  extensions  of  the  radial  field. 
More  precisely,  it  turns  out  that  in  many  cases  of  practical  interest,  the  above  expansion  converges  too  slowly 
to  be  competitive  with  the  resolution  of  the  full  system  of  equations,  which  we  discuss  now. 


Numerical  resolution  of  the  full  system  of  equations.  The  full  system  of  equations  correspond 
to  the  laws  of  conservation  of  mass  and  momentum.  In  the  simplest  case,  i.e.,  Mohr-Coulomb  materials, 
the  stress  equations  decouple  from  the  velocity  equations  (the  later  still  depend  on  the  stresses).  The  stress 
equations  correspond  then  to  a  pair  of  hyperbolic  conservation  laws  with  source  terms  and  nonlinear  boundary 
conditions.  To  date,  most  attempts  toward  the  numerical  simulation  of  the  corresponding  solutions  have  been 
based  on  a  judicious  change  of  variables.  Indeed,  the  yield  condition  corresponds  to  an  algebraic  constraint 
which  can  be  “solved”.  In  other  words,  when  using  the  so-called  Sokolowskii  variables  the  constraint  is 
automatically  satisfied.  Numerical  methods  can  be  then  used  to  solve  the  reformulated  system,  for  instance 
in  terms  in  the  Riemann  invariants,  which  was  our  first  approach. 

However,  and  in  spite  of  its  popularity,  the  use  of  those  Sokolowskii  comes  at  a  heavy  price;  the  con¬ 
servative  form  of  the  equations  is  lost,  and  thus,  so  is  the  ability  to  analyze  and/or  compute  shocks.  One 
has  to  conclude  that  the  equations  should  in  fact  preferably  be  solved  in  the  original  variables,  even  at  the 
price  of  having  to  deal  with  a  corresponding  complicated  hyperbolic  fiux  term  (no  a  priori  bounds  on  the 
wave  speeds).  Further,  the  discretization  of  the  equations  in  conservation  form  has  several  advantages  from 
a  purely  numerical  point  of  view.  We  have  implemented  an  algorithm  corresponding  to  a  formally  globally 
third  order  Discontinuous  Galerkin  method  for  the  resolution  of  three  dimensional  axisymmetric  hopper 
flows  involving  Coulomb  materials.  The  approach  has  been  extremely  successful. 


Computation  of  switch  stresses.  As  a  way  of  illustration,  we  include  hereafter  calculations  corre¬ 
sponding  to  the  influence  on  the  stress  field  of  abrupt  changes  in  the  wall  angle.  The  geometrical  situation 
is  as  follows 


Figure  1.  Geometrical  situation;  left:  transition  to  a  flatter  wall  angle,  right: 
transition  to  a  steeper  wall  angle.  The  radial  stress  is  used  to  generate  an  initial 
condition  on  the  curve  Fq.  The  domains  of  calculation  are  shaded. 

In  the  case  of  a  transition  from  a  steeper  to  a  flatter  hopper,  the  following  type  of  structure  is  observed 
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Figure  2.  Struciurt  of  the  stress  field  induced  by  a  transition  from  an  opening 
angle  o/lO"  to  one  0/15".  The  actual  geometry  of  the  hopper  is  represented. 

If  instead,  the  wall  goes  from  flatter  to  steeper,  the  stress  field  takes  the  following  form 


Figure  3.  Structure  of  the  stress  field  induced  by  a  transition  from  an  opening 
angle  o/lO"  to  one  0/12".  The  actual  geometry  of  the  hopper  is  represented. 

We  are  not  aware  of  other  successful  calculations  of  this  type  of  effects.  The  corresponding  equations 
are  quite  delicate  to  handle  numerically  because  of  the  interplay  between  the  source  terms  and  the  nonlinear 
boundary  conditions.  Further,  some  assumptions  related  to  the  structure  of  the  above  waves  can  be  found  in 
the  literature.  It  has  for  instance  been  conjectured  that  a  transition  to  a  flatter  hopper,  see  Figure  2,  would 
induce  shocks,  while  a  transition  to  a  steeper  one,  see  Figure  3,  would  corresponding  to  rarefaction  types 
of  waves.  Those  claims  are  typically  based  on  unsubstantiated  arguments,  for  instance  on  a  radial  solution 
analysis.  This  cannot  make  sense  here  because  of  the  lack  of  symmetry  of  the  full  domain,  see  Figure  1. 
Figures  2  and  3  suggest  that  the  situation  is  in  fact  more  complicated.  A  full  analysis  of  those  results  is 
under  way.  Experimental  comparisons  can  be  done  for  instance  based  on  the  values  of  the  normal  stress  on 
the  wall  predicted  by  the  numerics.  This  type  of  comparison  is  also  under  way,  specifically  for  cereal  types 
of  grains  in  steal  hoppers  for  instance. 
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Ongoing  and  future  work  about  granular  materials.  We  list  hereafter  ongoing  and  future  work, 
starting  from  tasks  that  are  well  under  way  or  expected  to  present  no  major  difficulties,  to  more  involved 
problems. 

(1)  Resolution  of  the  velocity  equations.  Within  our  current  model,  the  velocity  equations  consist 
in  an  additional  pair  of  hyperbolic  conservation  laws  coupled  to  the  stresses,  to  which  the  same 
numerical  approach  (essentially  the  same  Discontinuous  Galerkin  code)  can  be  applied.  One  of  the 
challenges  lies  here  in  determining  how  the  velocity  ^'information  travels” .  We  expect  a  careful 
numerical  study  to  bring  an  answer  to  this  type  of  questions,  that  are  for  now  open. 

(2)  Computation  of  flows  around  obstacles.  Assuming  the  flow  is  still  consistent  with  the  assump¬ 
tions  of  the  models,  this  computation  can  be  done  without  further  difficulties. 

(3)  Inclusion  in  the  models  of  interstitial  fluids.  When  dealing  for  instance  with  fine  powders  or 
soils  (clay,  etc...),  the  presence  of  interstitial  fluid  cannot  be  neglected.  The  models  have  then  to  be 
revised  to  include  the  corresponding  effects.  So  far,  our  work  in  this  direction  has  been  limited  to  the 
study  of  the  consolidation  of  fine  powders.  Further,  many  problems  of  prime  interest  to  the  Army 
such  as  mine  detection  or  soil  vehicle  interactions  are  by  nature  quite  close  to  both  the  computation 
of  Sows  around  obstacles  and  the  present  question  of  interstitial  Hows,  Plowing  is  another  example. 
Preliminary  contacts  have  been  taken  with  Drs.  J.F.  Peters  and  D.A.  Horner  at  the  U.S.  Army 
Engineer  Waterways  Experiment  Station  in  Vicksburg,  MS  (visit  to  WES  in  Oct.  98),  to  explore 
this  further. 

(4)  Determination  of  more  general  switch  stresses.  We  cannot  treat,  in  the  present  context,  the 
case  of  the  transition  from  a  vertical  pipe  to  a  conical  hopper.  The  reason  is  that  in  the  vertical 
pipe,  the  granular  material  is  thought  to  be  in  an  active  state,  as  opposed  to  a  passive  state  as  in  any 
conical  hopper.  This  change  corresponds  mathematically,  in  our  current  setting,  to  a  discontinuity 
in  the  hyperbolic  flux  function  itself.  The  model  as  it  stands  does  not  allow  for  the  prediction  of 
where  the  transition  between  active  and  passive  kind  of  flows  takes  place. 

(5)  Long  term  objectives.  Among  the  long  term  objectives,  we  include  the  possibility  of  solving 
the  above  problems  for  a  wider  range  of  constitutive  laws,  such  as  for  instance  critical  state  soil 
mechanics,  the  possibility  of  computing  and  analyzing  rigid-plastic  free  boundaries,  as  well  as  the 
treatment  of  fully  time  dependent  problems. 


4.  TECHNOLOGY  TRANSFER: 

We  will  pursue  our  collaboration  with  Jenike  &  Johanson,  Inc.  After  our  recent  visit  there,  the  focus  will 
be  on  the  efficient  resolution  of  the  full  problem  for  granular  flows  (as  described  above).  One  of  the  common 
goals  is  the  construction  of  numerical  solvers  that  work  in  nontrivial  geometries  (hopper  with  insert,  for 
instance),  and  for  a  wide  variety  of  materials.  Additional  contacts  have  already  been  scheduled.  Further, 
our  numerical  code  for  the  resolution  of  powder  consolidation  problems  in  currently  (Summer  99)  being  tested 
by  the  engineers  at  Jenike  &  Johanson.  This  code,  which  uses  a  Differential  Algebraic  Equation  approach 
to  the  problem  of  determining  stresses,  air  pressure  and  height  of  a  column  of  powder  in  an  axisymmetric 
geometry  is  of  immediate  practical  use. 


Figure  4.  Graphic  user  interface,  GUI,  for  the  powder  consolidation  problem, 
being  tested  at  Jenike  and  Johanson,  Inc. 

We  also  intent  to  explore  further  the  possibility  of  collaborating  with  the  researchers  at  the  U.S.  Army 
Engineer  Waterways  Experiment  Station  in  Vicksburg,  MS,  especially  in  terms  of  the  possible  use  the  present 
continuum  approach  to  “plowing  type  problems” . 
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